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PREFACE 


The  program  reported  herein  was  a joint  effort  between  Pacific- 
Sierra  Research  Corporation  (PSR)  and  Applied  Theory,  Inc.  (ATI).  This 
report  distills  only  the  positive  results  of  the  calculatlonal  program 
In  support  of  Test  Ban  Treaty  verification  and  does  not  discuss  the  many 
efforts  that  were  unfruitful.  The  authors  would  like  to  thank  their 
many  colleagues  at  PSR  and  ATI  who  contributed  to  the  program;  to 
Drs.  D.  Harkrider  and  E.  Herrin,  who  provided  valuable  insight  into 
the  problem;  and  to  Dr.  C.  Romney  and  Colonel  G.  Bui  in  at  ARPA/NMRO, 
who  helped  keep  us  on  track. 

This  Is  the  first  of  two  volumes  reporting  on  work  performed  under 
contract  to  the  Defense  Advanced  Research  Projects  Agency  from  October 
1975  through  November  1978.  This  volume  contains  the  executive  summary 
and  main  body  of  the  report.  Volume  II  contains  the  supporting  appendices. 
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EXECUTIVE  SUMMARY 


This  report  discusses  a calculational  program  aimed  at  improving 
U.S.  capability  to  verify  a Threshold  Test  Ban  Treaty  (TTBT)  by  seismic 
means . 

The  Soviets,  in  their  PNE  and  weapon  programs,  may  conduct  nu- 
clear tests  at  shallower  depths  than  the  United  States  or  the  French 
experience.  Since  yield  estimation  is  thus  required  at  depths  for 
which  little  or  no  direct  seismic  data  are  available,  we  emphasize 
shallow  bursts  in  our  analysis  and  examine  both  body-wave  and  surface- 
wave  effects. 

Two-dimensional  inelastic  source  calculations,  using  Applied 
Theory's  AFTON  program,  were  made  on  a representative  set  of  150-KT 
explosions.  This  set  covers  a range  of  depths  and  materials  that 
might  be  used  in  Soviet  weapon  and  PNE  programs.  These  sources  were 
used  to  develop  synthetic  seismometer  body-wave  signals  at  teleseismic 
distances.  Considerable  effort  was  devoted  to  excluding  nonpropagating 
waves  from  the  source  description. 

Simple  elastic  theory  indicates  that  the  primary  body  wave  from  a 
buried  explosion  receives  positive  or  negative  reinforcement  from  the 
free-surface  reflection--the  pP  wave.  The  amplitude  and  direction  of 
the  reinforcement  depend  on  timing,  and  therefore  on  depth  of  burial 
of  the  explosion.  Our  more  refined  calculations,  which  include  in- 
elastic and  gravitational  effects,  indicate  a "reflection"  amplitude 
smaller  than  that  calculated  for  the  elastic  case.  More  Important, 
the  reflected  wave  is  relatively  delayed  in  time,  so  that  transitions 
between  positive  and  negative  reinforcement  occur  at  shallower  depths 
of  bur1al--depths  usable  by  the  Soviets  and  therefore  of  crucial  im- 
portance in  verification  of  a 150-KT  TTBT.  Using  the  conventional  m^ 
definitions  and  the  AFTAC  Instrument  corrections,  we  find  the  depth 
of  burst  varies  with  mfe  by  at  least  ±0.2  in  the  150-KT  range. 

The  surface-wave  model  developed  is  based  on  Green's  function. 

Many  problems  were  encountered  In  modifying  the  AFTON  source-data  pro- 
gram to  provide  an  input  that  was  accurate  for  long-period  displacements- 
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and  to  extrapolate  data  well  beyond  the  reasonable  truncation  times 
for  the  program. 

Though  additional  calculations  must  still  be  performed,  our  in- 
vestigation has  yielded  the  following  preliminary  conclusions: 

1.  For  yield  determinations  and  comparisons,  inelastic  effects 
(including  cracking)  and  gravitational  effects  must  be  in- 
cluded in  the  analysis  of  either  body  waves  or  surface  waves 
from  shallow  explosions. 

2.  For  relatively  shallow  explosions  (say,  <1500  m at  150  KT), 
burial  depth  significantly  affects  both  the  body-wave  and 
surface-wave  signals.  This  effect  has  implications  for 
earthquake/explosion  discrimination  as  well  as  yield  estima- 
tion. 

3.  Without  correction  for  depth  effects  due  to  cracking  and 
gravity,  mb  (as  conventionally  defined)  can  result  in  sys- 
tematic yield  overestimation  of  shallow  Soviet  explosions  in 
the  150-KT  range  in  sedimentary  rocks. 

4.  For  burial  depths  between  200  m and  1500  m,  covering  the 
range  of  greatest  interest  for  TTBT  verification,  the  body- 
wave  yield  estimates  may  be  improved  by  using  uncorrected 
seismometer  amplitude  (A)  from  the  first  displacement  swing, 
or  from  the  largest  displacement  swing,  or  both--rather  than 
a combination  (conventionally  found  in  m^)  of  the  largest 
amplitude  with  its  associated  period  (T)  and  Instrument  cor- 
rection. Use  of  A,  alone  and  uncorrected,  also  appears  to 

make  m.  a more  direct  measure  of  radiated  energy.  Seismic 
b 

data  from  selected  explosions  should  be  reexamined  from  this 
viewpoint. 

5.  Topography  near  the  shot  point  does  not  appear  to  be  a major 
factor  in  seismic  signal  generation. 

Regarding  conclusion  4:  For  modern  seismometers,  output  ampli- 
tude is  nearly  proportional  to  the  speed  of  ground  motion  between  1 Hz 
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and  3 Hz,  the  frequencies  of  greatest  interest  for  explosion  monitor- 
ing. Thus,  the  approximate  logarithmic  proportionalities: 

log  (seismometer-amplitude)  ^ log  (ground-velocity)  % log  (ground-velocity )^. 

The  logarithm  of  the  largest  interval  spanned  by  peaks  and  troughs  of 

ground  velocity  is  roughly  proportional  to  the  logarithm  of  radiated 

kinetic  energy  (and  total  radiated  energy)--and  so,  therefore,  is  log  (A). 

On  the  other  hand,  at  the  frequencies  in  question,  A/T  represents  ground 

acceleration,  which  is  not  fundamentally  related  to  radiated  energy. 

A generalized  method  for  combining  body-wave,  surface-wave,  and 

other  measurements  to  refine  yield  estimates  has  been  reported  pre- 
★ 

vlously. 

Summaries  of  the  body-wave  and  surface-wave  calculations  to  date 
are  presented  in  the  body  of  this  report.  The  appendices  include  addi- 
tional details  and  numerical  methods. 

s 

A.  Ciervo,  et  al . , PSR  Notes  91,  93,  113,  and  122. 


C 


c 


c 

I 

I I 


t 


1 


preceding  page  blank  - not  filmed 


ix 

CONTENTS 

PREFACE  Hi 

EXECUTIVE  SUMMARY  v 

FIGURES  xl 

TABLES  xiii 

Section 

I.  THE  PROGRAM  OF  COMPUTATION  1 

Scope  of  the  Computational  Program  1 

Cratering  Calculations  and  Their  Significance  2 

The  Need  for  Calculation  to  Late  Problem-Times : 

Surface  Waves  4 

Special  Accuracy  Requirements  in  Computing  Body 

Waves  5 

Innovation  7 

II.  SURFACE  WAVES  9 

Introduction  9 

Source  Model  10 

Propagation  Model  13 

Results  16 

Conclusions  26 

III.  BODY  WAVES  27 

Body-Wave  Magnitudes  Versus  Shot  Depth:  Interaction 

of  Direct  and  Related  Waves  27 

Further  Comnents  on  Fig.  6:  Elastic  Versus 

Inelastic  Reflection  30 

Discontinuities  in  mb  35 

Specifics  of  Interference  and  Reinforcement: 

Approximation  of  Perfect  Reflection  36 

The  Factors  That  Enter  ma  and  mb,  and  Their 

Variation  With  Depth  42 

Correlation  of  Body-Wave  Source  Strength  With 
Features  of  Seismograms:  Further  Comments  44 

REFERENCES  47 


preceding  page  blank  - not  filmed 

xi 

FIGURES 


1.  Grid  Elements  Used  to  Model  Magic-Circle  Data  of  Source 

Volume  for  Wet  Sandstone  12 

2.  Schematic  Outline  of  PSR  Ms  Model  Showing  Source  Terms, 

Propagators,  and  Vertical  Response  for  One  Harmonic 
Component  and  One  Ring  Force  14 

3.  Surface-Wave  Magnitude  (M$)  Versus  Depth  of  Burial  in 

Granite  and  Wet  Sandstone  With  an  Average  Earth- 

Damping  Factor  Q = 300  17 

4.  Details  of  Shallowest  Granite  Simulation  22 

5.  Scatter  Plot  of  Ms  Versus  M2  for  All  Entries  in  Table  2 

(p.  18)  Showing  a Good  Correlation  25 

6.  Body-Wave  Magnitude  Versus  Depth  of  Burial,  Computed 

for  150  KT  Bursts  in  a Wet-Sandstone  Half-Space  28 

7.  Displacement  Pulses  in  the  Elastic  Near-Field  of  a 

150-Kiloton  Burst,  for  Three  Depths  of  Burial  in  a 
Half-Space  of  Wet  Sandstone  34 

8.  Ground  Motion  and  Seismometer  Response:  150-Kiloton 

Burst  in  a Wet-Sandstone  Half-Space;  Perfect  Ground- 
Surface  Reflection;  Three  Burial  Depths  38 

9.  Near-  and  Far-Field  Displacement  Pulses,  and  Seismometer 

Responses,  for  150-Kiloton  Bursts  at  Three  Depths  in  a 
Wet-Sandstone  Half-Space,  Calculated  With  Gravity  and 
Cracking  40 

10.  Near-  and  Far-Field  Displacement  Pulses,  and  Seismometer 

Responses,  for  150-Kiloton  Bursts  at  Three  Depths  in  a 

Wet-Sandstone  Half-Space,  Calculated  without  Gravity 

and  Cracking  41 

11.  Three  Depth-Dependent  Factors  That  Determine  Body-Wave 

Magnitudes  ma  and  mb.  Versus  Shot  Depth,  for  150- 

Kiloton  Bursts  in  a Wet-Sandstone  Half-Space  43 


A 


xiii 


preceding  page  blank  - not  filmed 


TABLES 

1.  Calculational  Matrix:  Tamped  150-KT  Explosions  3 

2.  Results  of  Nine  Simulated  Explosions  With  Sensitivity 

Test  18 

3.  Amplitude  and  Duration  of  the  Relief  Wave  Due  to  a 

150-KT  Burst  in  a Wet-Sandstone  Half-Space  Versus 

Depth  of  Burial  32 


I.  THE  PROGRAM  OF  COMPUTATION 


The  main  goal  of  the  program  reported  herein  was  to  determine,  more 
accurately  than  before,  how  teleseismic  surface  and  body  waves  driven  by 
buried  explosions  vary  with  shot  depth.  At  first,  emphasis  fell  on  shallow 
Soviet  PNEs,  and  the  scope  of  effort  later  widened. 

A computer  code  (AFTON  2A)  for  calculating  fields  of  axisymmetric 
continuum  motion  was  the  key  tool  used  to  define  the  ground  motions 
studied.  Since  actual  depths  of  burial  seldom  exceed  a kilometer  (km), 
and  many  buried  bursts  have  produced  craters,  it  was  necessary  in  our 
computational  work  to  take  account  of  tensile  cracking  and  gravity. 

SCOPE  OF  THE  COMPUTATIONAL  PROGRAM 

At  the  outset,  we  hoped  to  compute  fields  of  motion  for  18  buried 
bursts.  Early  on,  that  number  was  reduced  to  12  when  it  was  found  that 
inelastic  processes  (mainly  cracking  with  gravity)  lengthened  by  at  least 
a factor  of  two  the  time  of  significant  motion  in  the  near  field  (i.e., 
within  a few  kilometers  of  the  shot  point  for  bursts  yielding  150  kilotons 
(KT)  of  energy).  The  hypothetical  host  media  for  the  12  bursts--each  of 
150-KT  yield--were  Piledrlver  granite  (3  shot  depths),  wet  sandstone  (5 
depths),  dry  sandstone  (3  depths),  and  weak  dry  sandstone  (1  depth). 

Granite  and  sandstone  were  chosen  as  types  of  material  in  which  Soviet 
PNEs  would  be  likely  to  occur.  Degree  of  saturation,  a key  property  of 
earth  media  that  affects  both  strength  and  compactibility,  was  varied 
(from  0%  to  100%)  in  the  change  from  wet  sandstone  to  dry.  Shear  strength, 
another  key  property,  was  also  varied;  its  effects  on  ground  motion  due 
to  PNEs  were  explored  via  the  change  from  dry  sandstone  to  weak  dry  sand- 
stone. Similarly,  the  effects  of  a variation  in  dry  porosity  alone  (from 
10%  to  0%)  are  seen  in  the  fields  computed  for  weak  dry  sandstone  and  wet 
sandstone,  whose  strengths  were  identical. 

The  12  bursts  covered  5 different  depths,  viz.,  (10,  30,  39,  47.622, 
1/3 

and  100)  x W ' meters  (m),  where  W is  the  yield  in  kilotons.  Containment 
was  expected  only  for  the  depth  100  Wly/3,  and  shots  at  depths  of  10  W1//3 
and  100  W ' m were  simulated  only  in  wet  sandstone.  Also,  since 
30  x 4^3  = 47.622,  the  150-KT  field  for  depth  30  x 150^3  m is  a given 
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medium  provided  a simply  scaled  version  of  ground  motion  from  a 600-KT  burst 

1/1 

at  depth  47.622  x 150  m.  That  600-KT  field  could  be  compared  with  the 
150-KT  field  computed  for  the  same  depth  (47.622  x 1501/3  m).  A similar 
comparison  was  possible  between  37.5-KT  and  150-KT  shots  at  a depth  of 
30  x 150  ' m.  In  addition,  to  assess  gravity-caused  departures  from 
simple  scaling  rules,  a calculation  was  made  of  motion  for  a yield  of 
600  KT  at  a depth  of  30  W^3  m.  Then,  as  interest  in  Soviet  test-media 
grew,  results  obtained  in  another  program  for  shots  in  a "Soviet  granite" 
were  augmented  by  three  further  calculations  for  that  medium--two  at  depths 
other  than  the  five  already  noted.  Most  important,  at  the  very  end  of  the 
program,  motion  was  computed  for  a 150-KT  shot  at  a depth  of  744  rr,  in  wet 
sandstone.  (Table  1 summarizes  the  full  problem  set.) 

The  last  calculation  cited  was  one  of  a series  of  nine,  made  necessary 
bv  results  that  surfaced  at  a late  stage  of  the  program.  The  nine  final 
calculations  were  made  after  modifying  the  AFT0N-2A  code  somewhat,  in  ways 
that  improved  its  accuracy  and  thereby  made  it  possible  to  determine  ground 
motion  at  greater  depths  and  later  times  than  before.  Besides  the  shot 
744  m deep  in  wet  sandstone,  the  final  series  included  calculations  of 
motion  for  the  three  granite  bursts  of  Table  1,  and  the  two  in  wet  sand- 
stone at  depths  of  39  W^3  and  100  W1,/3  (207.2  m and  531  m)--all  five  of 
which  had  been  run  earlier.  Moreover,  each  of  the  three  wet-sandstone 
problems  of  the  final  series  was  run  twice:  once  with  gravity  and  crack- 
ing; once  without. 

The  motivation  for  the  final  nine-problem  series,  and  the  code  modi- 
fications that  preceded  it,  are  discussed  below.  Here,  we  note  that 
relative  values  of  body-wave  magnitude  (m^)  originally  obtained  for  granite 
changed  by  no  more  than  0.03;  for  the  two  recomputed  wet-sandstone  bursts, 
relative  m^-values  changed  by  0.22--but  only  because  of  slight  differences 
In  seismometer  response  that  bridged  an  mb-discontinuity  (see  Sec.  III). 

CRATERING  CALCULATIONS  AND  THEIR  SIGNIFICANCE 

In  the  past,  the  dimensions  of  craters  calculated  for  nuclear  surface 
bursts  using  codes  like  AFT0N  2A  have  been  smaller  by  factors  cf  10  to  50 
than  observations  of  actual  craters  suggest.  Yet,  it  is  also  true  that 
the  codes  predict  realistic  craters  for  buried  explosions  [1]  even  when 
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shot  points  are  much  shallower  than  the  depths  of  interest  herein  ( >50  m 
for  150  KT).  Those  results  are  explained,  at  least  in  part,  by  the 
simplicity  of  the  process  of  coupling  explosive  energy  to  the  ground,  when 
shot  depths  exceed  ^3  m (for  150  KT).  Also,  unless  shot  depths  are  large 
(on  the  order  of  1 km  for  a 150-KT  burst),  elastic  tensile  stresses  above 
the  shot  point  so  greatly  exceed  strength  limits  that,  to  predict  cracking, 
little  detailed  knowledge  is  needed  either  of  tensile  strengths  or  of 
cracking  mechanisms.  Accordingly,  as  an  important  result  of  our  cratering- 
burst  calculations,  we  were  able  to  show  that  an  adequate  account  of  the 
effects  of  cracking  and  gravity  is  given  by  AFTON  2A.  Specifically,  the 
volume  of  our  calculated  Sedan  Crater  was  60-percent  greater  than  that 
observed  (versus  error-factors  of  10  to  50  for  surface  bursts);  and,  for 
granite  and  wet  sandstone,  our  calculated  dependence  of  crater  size  on 
shot  depth  was  consistent  with  observation. 

The  ability  to  compute  crater  sizes  was  comforting,  but  insufficient, 
to  secure  the  main  goal  of  the  program.  Rather,  our  toughest  computational 
problems  arose  from  a need  to  define  fields  of  motion  at  much  later  times 
(>2  sec  for  150-KT  bursts)  than  those  to  which  explosion  calculations  are 
usually  run;  by  contrast,  for  buried  explosions,  crater  dimensions  are 
largely  determined  by  the  early  history  of  ground  velocity  (<0.5  sec  for 
150-KT  bursts  in  granite  and  wet  sandstone). 

THE  NEED  FOR  CALCULATION  TO  LATE  PROBLEM-TIMES:  SURFACE  WAVES 

At  about  the  time  that  the  last  of  the  program's  three  one-year  phases 
began,  a major  result  emerged:  Surface-wave  magnitudes  were  found  to  be 
roughly  proportional  to  the  logarithm  of  the  net  downward  impulse  delivered 
by  an  explosive  source  to  the  earth  (Sec.  II).  More  precisely,  the  critical 
quantity  proved  to  be  the  vertical  impulse  delivered  to  the  earth,  in  a few 
tens  of  seconds,  across  any  surface  S in  the  elastic  near-field  of  the 
burst--provided  that  the  shot  point  is  enclosed  by  S-pl us-the-ground-sur- 
face.  The  surface  S was  taken  as  a hemisphere  H,  viz.,  the  lower  half  of 
a sphere  with  its  center  at  ground  zero,  and  having  about  the  smallest 
radius  consistent  with  elastic  deformation  (^1.2  km  for  wet  sandstone  and 
2. 1 km  for  granite) . 

Apart  from  in-flight  ejecta  for  the  shallowest  shot  depths  studied, 
all  significant  motion  of  material  originally  within  the  surface  H ended 
before  10  sec--from  which  it  follows  that  the  vertical  impulse  delivered 
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across  H in  10  sec  is  zero.  Unfortunately,  in  all  finite-difference  or 
finite-element  codes,  the  use  of  discrete  analogs  of  the  exact  equations 
of  motion  to  effect  their  numerical  integration  is  attended  by  discretiza- 
tion error.  As  a result,  the  vertical  impulse  found  from  our  first  calcula 
tions  was  not  exactly  zero;  that  spurious  impulse,  though  small  in  terms 
of  the  discretization  error  required  to  produce  it,  dominated  our  first 
values  of  surface-wave  magnitude,  M$. 

For  our  purposes,  the  net  effect  of  discretization  error  was  to  im- 
part a nonphysical  acceleration  of  ^0.002  g to  ground  material  in  our 
"static  equilibrium"  fields  (g  = gravitational  acceleration).  In  most 
applications,  an  error  of  that  size  would  be  negligible;  in  the  present 
context  it  was  unacceptably  large.  Translated  into  vertical  displacement, 
for  example,  it  amounted  to  about  6 centimeters  (cm)  after  2.5  sec  of 
motion  (the  final  time  of  most  of  our  calculations);  by  comparison,  the 
peak  displacement  for  wet  sandstone  was  ^24  cm  at  the  bottom  of  the  hemi- 
sphere H.  Thus,  the  error  was  big  enough  to  affect  slightly  the  first 
values  computed  for  mb,  although  it  was  serious  only  for  M$.  As  discussed 
below,  however,  it  became  critical  for  mb  as  well,  when,  late  in  the 
program,  new  facts  came  to  light  as  to  the  properties  of  elastic  near- 
fields. In  the  end,  to  compute  surface-  and  body-wave  fields  reliably 
enough  for  this  program,  the  simplest  course  was  to  eliminate  the  error 
entirely  (see  pp.  7-8).  That  was  done,  and  5 of  our  original  12  problems 
were  then  rerun;  a new  one  was  also  run. 

SPECIAL  ACCURACY  REQUIREMENTS  IN  COMPUTING  BODY  WAVES 

A major  result  of  our  fif*st  calculations  above  was  the  discovery  that 
ground  displacements,  as  seen  near  the  bottom  of  the  hemisphere,  H, 
varied  greatly  with  depth  of  burial.  They  were  also  much  different  in  thei 
late  stages  (after  ^1  sec  of  motion)  than  we  were  led  to  expect  on  the 
basis  of  (inexact)  elastic  half-space  theory.  Further,  the  duration  and 
timing  of  surprising  displacement-pulse  features  were  such  that  tensile 
cracking  (which  was  widespread,  according  to  the  calculations),  and 
gravitational  subsidence  of  the  cracked  medium,  could  have  produced  them. 

In  late  1977,  the  approximate  picture  of  elastic  near-field  motion 
previously  available  gave  way  to  a more  accurate  one,  whose  development 


6 


must  be  counted  a major  advance  in  the  understanding  of  elastic  near-fields 
(surprisingly)  [2]:  The  exact  near-field  of  a compressive,  burst-like 
point-source  in  an  elastic  half-space  was  derived  (or  rederived  [3]  ). 

Thus,  early  in  1978,  we  learned  that  late-time,  near-field  displacements 
similar  to  those  seen  in  our  calculations  would  occur  even  in  a gravity- 
free,  uncrackable  elastic  medium.  Furthermore,  those  elastic  displace- 
ments would  decay  so  rapidly  with  distance  from  the  shot  point  as  not  to 
affect  teleseismic  body  waves;  they  are  "local  parts"  of  the  complete 
near-field.  More  precisely,  drawing  on  the  principle  of  superposition, 
the  parts  of  a linear  near-field  whose  absence  would  leave  the  far  field 
unaffected,  are  referred  to  below  as  "local";  the  part  of  the  near  field 
that  cannot  be  altered  without  changing  the  far  field  is  called  "global." 

The  new  knowledge  of  elastic  near-fields  plainly  made  it  necessary 
to  distinguish  between  elastic  contributions  to  our  calculated  near-field 
displacements  and  all  other  contributions.  In  short,  our  key  task  became 
that  of  determining,  as  a function  of  burial  depth,  the  influence  on 
teleseismic  body  waves  of  processes  not  amenable  to  simple  theoretical 
treatment--particularly , cracking  and  gravitational  subsidence. 

To  separate  elastic  displacements  from  all  others  in  the  near-field 
of  a given  burst,  we  found  it  easiest  to  perform  two  calculations  that 
were  identical,  except  for  the  inclusion  of  gravitational  acceleration  and 
cracking  in  one  of  them,  and  not  in  the  other.  In  both  instances  (since 
body-wave  fields  are  irrotational ) , only  the  irrotational  part  of  each 
near-field  displacement  pulse  was  wanted.  That  requirement,  given  avail- 
able computational  tools,  was  best  met  by  calculating  motion  at  much 

3r 

The  effect  of  local  terms  on  elastic  near-fields  was  recognized 
first  by  S.  Day  (T.  Bache;  telephone  discussion  in  October  or  November 
1977).  A check  of  the  perfect-reflection  approximation  (Sec.  Ill) 
against  an  exact  solution  was  called  for  on  7 March  and  8 August  1977 
by  J.  Trulio  (at  NMRO-sponsored  meetings),  who  suggest  d that  an  S3 
elastic-wave-propagation  code  be  used  to  get  it  (analytic  methods  were 
used  instead). 

★ ★ 

Reference  3 is  based  on  a solution  to  the  problem  of  a compressive 
(spherical)  time-varying  source  in  a half-space  due  to  L.  Cagniard,  in 
the  translation  of  Flynn  and  Dix  ( Reflection  and  Refraction  of  Progressive 
Seismic  Waves,  McGraw-Hill,  1962)  as  corrected  by  J.  Workman.  Although 
surface  waves  are  omitted,  the  results  in  AFSWC-TR-68-27  are  virtually 
Indistinguishable  from  the  complete  exact  solution  at  steep  takeoff  angles, 
at  least  on  the  hemisphere  H (and  below). 
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greater  depths  than  in  the  past.  Shear  and  surface  waves  were  thereby 
removed  from  the  displacement  pulses  that  define  explosions  as  sources 
of  teleseismic  signals,  simply  because  these  waves  travel  too  slowly  to 
contribute  to  the  first  few  seconds  of  deep  near-field  motion. 

To  obtain  displacement  pulses  deep  in  the  ground,  motion  had  to  be 
computed  to  even  later  problem-times  than  before  (>3  sec,  it  was  found). 
Hence,  as  in  characterizing  buried  bursts  as  sources  of  teleseismic  sur- 
face waves,  forces  again  had  to  be  in  exact  equilibrium  in  the  static 
AFT0N-2A  field.  For  example,  at  the  relatively  great  depths  reached  in 
our  final  calculations  (^6  km),  peak  displacements  were  enly  ^2  cm, 
whereas  a discretization  error  of  0.002  g In  acceleration  would  have 
altered  late-time  displacements  by  10  cm. 

A computer  code  was  available  during  the  program  that  was  supposed 
to  propagate  body  waves  correctly  to  teleseismic  distances,  starting 
from  the  hemispherical  near-field  surface  H,  even  if  the  global  part  of 
the  irrotational  field  came  mixed  with  its  local  parts,  and  with  shear 
and  surface  waves.  However,  the  effects  of  cracking  and  gravity  on  near- 
field motion,  free  of  shear  and  surface  waves,  and  largely  free  of  local 
contributions,  could  only  be  determined  by  calculating  the  near-field 
at  great  depths  and  late  times.  To  do  so,  it  was  necessary  to  push  the 
computational  art  forward  a bit,  as  described  below.  When  that  was  done, 
we  found  that  the  large  effects  on  body  waves  attributed  earlier  to  gravity 
and  cracking,  both  in  the  near  field  of  a buried  burst  teleseismical ly, 
remained  large  (Sec.  III). 

INNOVATION 

The  unique  flexibility  of  the  AFT0N-2A  code  makes  it  practical  to 
compute  motion  from  a 150-KT  explosion,  without  assumptions  of  small  dis- 
placement, linearity,  etc.,  and  to  depths  as  great  as  6 km  and  times 
longer  than  3 sec.  Specifically,  a domain  of  AFT0N-2A  computation  can  be 
held  fixed  in  size  while  being  moved  continually  downward  at  the  speed  of 
longitudinal  waves.  The  cost  of  advancing  a numerical  solution  through 
unit  time  is  then  constant,  and  the  total  cost  of  computation  is  pro- 
portional to  the  time  after  burst  at  which  calculation  is  stopped.  On 
the  other  hand,  if  all  disturbed  material  had  to  be  included  in  the 
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computational  domain  (usual  with  computer  codes),  then  the  cost  of  calcula 
tlon  would  increase  roughly  as  the  cube  of  a final  time  of  calculation, 
rather  than  the  first  power. 

To  achieve  exact  static  equilibrium  in  a discrete  field  proved  more 
difficult  than  to  avoid  excessive  computation.  As  part  of  the  overall 
discretization  process,  a continuous  medium  is  divided  into  finite  poly- 
hedral cells  of  material,  whose  vertices  form  a "finite  difference  mesh" 
[4].  By  adopting  a simple  approximation  to  the  structure  of  an  exact 
continuum  field  within  a cell,  discrete  equations  of  motions  are  produced. 
The  properties  of  those  equations--and  particularly  the  absence  of  pro- 
perties not  possible  for  a physical  conti nuum--depend  on  the  local  (cell- 
scale)  field-structure  assumed.  The  properties  of  the  discrete  equations 
can  be  made  completely  consistent  with  those  of  the  continuum  equations 
they  replace,  but  only  if  stress  and  strain  are  locally  homogeneous  (for 
reasons  that,  though  not  complicated,  would  take  us  too  far  afield  to 
present  here). 

Some  of  the  unique  and  basic  transformation  properties  of  the  AFTON 
equations  of  motion- -which  include  exact  conservation  of  mass,  momentum, 
and  energy--are  tied  to  the  assumption  of  homogeneous  strain  and  stress 
fields  within  cells.  On  the  other  hand,  overburden  stresses  in  a static 
equilibrium  field  (under  gravity)  increase  continuously  with  depth.  Those 
two  different  field-structures  in  the  AFTON  equations  were  reconciled-- 
but  the  effort  was  not  trivial.  Specifically,  while  preserving  all 
transformation  properties  of  the  AFTON  equations,  as  well  as  the  nearly 
complete  freedom  of  spatial  mesh-point  distribution  offered  by  AFTON  2A, 
the  code  has  now  been  modified  so  that  its  gravitational  and  material 
forces  cancel  exactly  in  static-equilibrium  fields.  So  modified,  the 
code  was  used  to  recompute  the  three  granite  fields  studied  in  the 
program,  and  two  of  the  original  wet-sandstone  fields  (pp.  1-2,  above). 
Also  for  wet  sandstone,  motion  was  computed  for  one  new  shot  depth,  and 
fields  were  calculated  for  the  last  three  explosive-sandstone  systems 
both  with  and  without  cracking  and  gravity. 


II.  SURFACE  WAVES 


INTRODUCTION 

The  PSR  M$  model  is  based  on  Harkrider's  analytic  expressions  for 
surface-wave  Green's  functions  in  a flat  multilayered  half-space  [5].  The 
multilayer  model  of  the  earth  structure  is  adopted  because  the  modelled 
explosion  sources  with  which  we  are  concerned  here  all  have  cylindrical 

symmetry,  and  the  symmetry  of  the  layered  earth  model  is  also  cylindrical. 

3 

Previous  efforts  by  S (Systems,  Science  and  Software  [6])  to  calculate 
M$  from  the  same  source  model  failed  partly  because  of  the  difficulty  in 
finding  an  equivalent  point-source  term  for  use  with  the  whole-space  Green's 
function.  In  fact,  there  is  no  equivalent  point-source  term,  because  the 
entire  free  surface  becomes  a source  when  whole-space  Green's  functions 
are  used.  The  PSR  M$  model  uses  half-space  Green's  functions  and  there- 
fore avoids  the  source  terms  from  the  free  surface.  This  method  results 
In  a source  that  is  localized  near  the  explosion  and  therefore  easier  to 
treat  mathematically. 

The  following  equation  used  to  calculate  M$,  derived  from  Marshall  and 
Basham  [7],  is  dependent  only  on  the  maximum  vertical  ground  displacement 
near  the  20-sec  periods  at  a distance  of  10,000  km: 

M$  = 9.2  + log  ( Umax [ cm] ) . (1) 

The  constant  9.2  has  the  distance  and  period  correction  already  included 
according  to  the  prescription  in  Ref. '7.  The  value  of  Umax  is  taken  to  be 
the  maximum  displacement  from  equilibrium  in  either  the  positive  or 
negative  direction. 

The  following  discussion  of  the  PSR  M$  model  is  directed  toward  a 
complete  specification  of  the  source  parameters  used,  followed  by  a 
description  of  the  propagation  method  of  calculating  Um_  for  sources  of 
Interest.  Note  that  this  model  calculation  makes  several  simplifying 
physical  assumptions,  which  are  discussed  briefly  as  necessary. 
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SOURCE  MODEL 

The  following  details  of  the  nonelastic  source  model  constitute  a 
mathematical  description  of  the  source  suitable  for  computing  surface-wave 
magnitudes.  Basically,  we  are  dealing  with  a two-dimensional  problem, 
since  the  theoretical  problem  of  a spherical  explosion  in  a homogeneous 
isotropic  half-space  has  cylindrical  symmetry.  The  plastic-flow  character  j 

of  the  source  need  never  enter  into  the  description  of  the  source  vis-a-vis 
surface-wave  creation.  The  volume  of  space  where  plastic  flow  and  all 
other  nonelastic  behavior  takes  place  can  be  surrounded  by  an  imaginary 
surface  outside  of  which  the  usual  equations  for  small -amplitude  linear 
elastic-wave  propagation  are  valid.  By  surrounding  the  source,  only  the 
specification  of  the  dynamic  state  of  stress  for  material  in  contact  with 
the  imaginary  surface  is  relevant  for  the  propagation  calculation,  and 
the  details  of  rock  mechanics,  etc.,  which  produced  that  dynamic  state, 
are  no  longer  needed  in  the  source  description.  Of  course,  by  following 
this  source  description,  we  could  have  sources  that  are  not  physically 
possible.  The  above  definition  of  any  particular  source,  however,  will 
have  certain  average  values  that  give  clues  to  its  faithfulness  to  physical 
reality;  e.g.,  net  vertical  impulse  should  be  zero.  In  practice,  the 
imaginary  surface  on  which  the  elastic  source  is  specified  has  been  a 
closed  numerical  grid  of  vertical  cylindrical  strips  and  horizontal  annular 
strips,  located  about  2-km  away  from  the  shot  point.  At  one  point  during 
this  study,  the  imaginary  surface  was  a circular  arc  that  was  rotated  about 
the  axis  of  symmetry  to  produce  a hemisphere  centered  at  ground  zero.  This 
hemisphere  acquired  the  name  "magic  circle,"  which  is  still  applied  to  the 
numerical  grid  although  it  is  no  longer  strictly  circular. 

The  complete  mathematical  specification  of  magic-circle  data  needed 
for  a formal  propagation  calculation  can  be  determined  by  looking  at  the 
wave  equation  as  an  initial  value  problem  in  the  theory  of  partial  differ- 
ential equations.  By  analogy  with  the  one-dimensional  wave  equation,  the 
value  of  the  scalar  potential  and  its  first  derivative  must  be  specified 
on  the  boundary  point  to  completely  determine  the  solution  everywhere. 

The  analogous  boundary  conditions  for  the  source  considered  here  are  the 
values  and  derivatives  of  the  potential  over  a two-dimensional  surface 
embedded  in  three  space  dimensions  and  one  time  dimension.  In  practice, 
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we  would  have  to  specify  the  potential  and  its  derivatives  over  the  hemi- 
sphere for  a long  time  duration.  The  PSR  M$  model,  however,  is  limited 
1)  because  it  uses  force,  not  potentials;  and  2)  because  derivatives  are 
not  used  at  all. 

The  source  description  currently  in  use  is  a collection  of  directed 
ring  forces  located  on  the  hemisphere  swept  out  by  the  magic  circle  and 
with  time  dependencies  determined  by  the  numerical  results  of  the  modelled 
explosions  as  calculated  by  Applied  Theory,  Inc.  We  adopt  the  simplified 
source  model  for  two  reasons:  1)  the  solution  is  much  simpler  to  implement, 
since  point-force  Green's  functions  were  found  in  the  literature;  2)  good 
test  cases  have  not  yet  been  developed  for  the  purpose  of  checking  even 
the  simple  solution,  and  so  the  more  complicated  one  would  be  premature 
until  such  a testing  ability  is  developed,  since  the  simpler  method  may  be 
verifiable  for  explosion  sources. 

It  is  possible  to  get  some  idea  of  the  validity  of  using  a collection 
of  ring  forces  by  trying  to  imagine  what  the  actual  elastic  solution  for 
the  ring  forces  would  be,  compared  to  the  continuous  distributed  force 
model.  First,  the  singularity  in  stress  near  the  ring  force  is  nonphysical. 
But  If  one  were  to  observe  the  motion  at  a distance  from  the  rings  about 
equal  to  the  ring  spacing,  then  it  is  easily  imagined  that  the  ring-force 
solution  would  begin  to  resemble  the  distributed-force  solution,  although 
the  details  would  depend  on  interference  effects  between  the  rings.  In 
particular,  it  seems  that  ring  forces  would  radiate  energy  both  inward  and 
outward  equally,  at  least  locally,  whereas  distributed  forces  would  radiate 
energy  predominantly  outward  relative  to  the  magic  circle.  Intuitively, 
one  would  expect  that  if  the  response  due  to  the  ring  forces  at  a moderate 
distance  outside  the  magic  circle  could  be  made  by  diffraction  to  resemble 
the  distributed-force  response,  then  the  same  diffractive  effect  would 
limit  energy  radiation  inward  as  well.  But  even  if  diffraction  effects  did 
not  operate  as  hoped,  the  energy  budget  would  still  not  be  off  by  a factor 
of  more  than  two;  and  as  long  as  the  different  explosions  do  not  differ 
radically,  which  is  the  case,  the  relative  values  of  Ms  should  still  be 
meaningful . 

Figure  1 represents  the  surface  elements  for  a wet-sandstone  grid  on 
which  the  ring  forces  are  defined.  Note  that  the  surface  elements  are  all 


components  are  shown  for  one  of  the  rings  at  one  instant  of  time.  DOB  = depth  of  burial. 
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nearly  vertical  or  nearly  horizontal,  but  the  direction  of  force  on  either 
does  not  have  to  align  with  the  surface  normal.  There  are  typically  about 
60  directed  force  elements  for  each  explosion  event,  and  the  time  duration 
runs  about  2 sec. 

PROPAGATION  MODEL 

The  basic  element  of  the  propagation  model  is  the  surface-wave  Green's 
function  that  describes  in  one  complex  number  the  teleseismic  vertical 
ground  motion  (amplitude  and  phase)  when  the  source  is  a buried  point  force 
(horizontal ly  or  vertically  directed)  that  has  harmonic  time  dependence. 
Figure  2 represents  what  the  Green's  function  does--which  is  to  solve  for 
the  teleseismic  vertical  displacement  due  to  one  component  of  our  model 
source.  The  solution  for  the  complete  source  is  found  by  superposition 
in  time  and  space  of  the  component  parts. 

The  desired  result  is  the  vertical  ground  motion  U(t)  (in  centimeters) 
at  a large  but  arbitrary  range  of  10,000  km  as  a function  of  time.  Once 
this  function  is  found,  the  maximum  value,  Umax,  is  used  in  Eq.  (1)  to  give 
M$.  The  full-time  dependence  is  found  in  the  usual  way  by  Fourier  synthesis 
of  all  component  harmonic  terms: 


U(t)  = 


co 


— oo 


(2) 


The  Fourier  transformed  displacement,  U(cd),  depends  on  the  source 
terms  only  at  the  angular  frequency  id;  but  on  the  individual  force  elements 
located  at  all  different  positions  on  the  magic  circle.  The  teleseismic 
response  U j (u>)  for  the  j-th  ring  force  is  a linear  combination  of  both 
the  horizontal  and  vertical  components  of  the  force  element  at  the  fre- 
quency id,  Pyj(u>),  and  Pj_j j (id)  , multiplied  by  the  horizontal  and  vertical 
ring-force  Green's  functions,  GVjU)  and  (u>) : 


Uj(u.)  = PyjGyj  + pHjGHj 


(3) 
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The  force  terms  and  P^.  are  the  Fourier  components  of  each  ring 
force  at  the  angular  frequency  w.  In  general,  the  time  duration  of  the 
force  for  the  simulated  explosion  events  runs  out  to  about  2 sec,  whereas 
the  frequencies  of  interest  correspond  to  about  20-sec  periods.  Since  any 
finite  pulse  has  all  frequency  components  included,  the  individual  2-sec 
force  pulses  will  have  20-sec  components  as  well,  but  of  very  small 
amplitude.  At  first,  one  is  tempted  to  assume  that  e1a)t  changes  so  little 
for  0 < t < 2 sec  for  u>  corresponding  to  20-sec  periods  that  the  20-sec 
component  is  nearly  the  same  as  the  "zero-frequency"  component.  In  general, 
that  would  be  a good  approximation.  But  for  the  particular  mechanical 
situation  here,  the  conservation-of-momentum  principle  requires  that  the 
net  vertical  impulse  transmitted  by  the  magic  circle  should  be  equal  to 
zero.  Therefore  the  sum  over  all  vertical  force  elements  has  a zero 
frequency  component  equal  to  zero: 


(4) 


In  practice,  Eq.  (4)  is  beyond  the  range  of  the  modelled  sources  be- 
cause the  lowest  frequencies  are  not  strictly  zero,  but  instead  correspond 
to  2-sec  time  periods;  and  at  those  periods,  the  conservation  of  momentum 
applied  at  the  magic  circle  is  affected  by  ejecta  that  may  still  be  in 
freefall  for  some  explosion  events. 

We  therefore  must  be  somewhat  cautious  in  using  the  small -phase 
approximation  for  e”1ut,  since  we  know  the  20-sec  component  and  the  zero- 
frequency  component  of  the  total  vertical  impulse  are  not  equal.  In  the 
propagation  calculation,  we  have  adopted  the  small-phase  approximation, 
and  have  forced  Eq.  (4)  to  be  satisfied  by  modifying  the  magic-circle 
data  in  a simple  manner  that  simulated  the  effect  of  all  returning  ejecta. 

The  last  piece  to  the  outline  of  the  propagation  calculation  is  the 
pair  of  ring-force  Green's  functions  Gyj  and  . The  force  functions 
are  calculated  from  the  point-force  Green's  functions  g^  and  g^  of  Ref.  5 
for  the  harmonic-force  components,  corresponding  to  the  depth  and  radius 
of  the  j-th  ring  force  of  the  magic  circle.  The  ring  functions  Gy  and  G^ 
are  found  from  the  point  functions  gy  and  gH  by  integration  of  the  point- 
force  density  per  radian  over  a complete  circle  of  the  azimuth  angle. 
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Note  that  all  contributions  of  the  azimuthal  integration  are  additive  for 
the  vertical  force,  but  there  is  cancellation  on  opposite  sides  of  the 
circle  for  the  horizontal  force.  This  cancellation  results  in  a ratio 
G^/Gy,  which  is  reduced  from  the  ratio  q^/g^  by  a factor  of  about  nr/A, 
where  r is  the  radius  of  the  ring  (*sl  km)  and  A is  the  Rayleigh  wave- 
length («60  km). 

The  point-force  Green's  functions  used  here  differ  from  those  of 
Harkrider  [5]  in  two  relatively  minor  ways.  First,  the  actual  multilayered 
Green's  functions  are  replaced  by  an  approximation  in  which  both  the  ampli- 
tude and  the  depth  dependence  are  adjusted  for  each  frequency  so  that  it 
looks  like  a homogeneous  half-space  Rayleigh  wave.  This  adjustment  makes 
use  of  the  dispersion  relation  for  the  multilayered  half-space,  but  it 
avoids  the  necessity  of  doing  the  full  multilayered  calculation  once  the 
dispersion  relation  is  known.  The  dispersion  relation  published  by 
McEvilly  [3]  for  North  America  is  used  here.  Second,  an  amplitude-damping 
model  is  added  so  that  a better  estimate  of  absolute  value  of  M can  be 
obtained,  although  it  is  still  felt  that  relative  values  of  M$  are  more 
meaningful.  The  damping  model  uses  a constant  Q damping  factor,  where  Q 
is  defined  as  the  inverse  of  the  fractional  energy  lost  per  cycle  of  a 
wavetrain.  The  justification  for  a constant  Q model  is  given  by  Knopoff 

[9] ;  and  the  value  Q = 300,  which  is  used  in  this  model,  represents  a 
compromise  between  values  suggested  by  Knopoff  (Q  = 150)  and  Tryggvason 

[10]  (Q  = 1300). 

RESULTS 

Results  presented  in  this  section  are  of  three  basic  types:  1)  the 
versus  depth  results  of  the  physical  model;  2)  a sensitivity  study  of 
the  effects  of  changing  a few  numerical  parameters  on  the  last  data  frame 
of  the  time  sequence;  3)  an  identification  of  relevant  source  parameters 
that  correlate  strongly  with  M . 

Figure  3 and  Table  2 show  the  M$  versus  depth  results  for  the  nine 
explosion  events  studied.  Of  the  nine,  all  three  shots  in  granite  re- 
sulted in  craters,  whereas  only  the  shallowest  one  in  wet  sandstone 
cratered.  The  shallowest  one  in  "elastic"  wet  sandstone  would  have 
cratered  had  not  the  material  been  artifically  prevented  from  failing  in 
tension.  Table  2 shows  a list  of  results,  including  cratering,  for  the 
shots  studied. 
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Table  2 

RESULTS  OF  NINE  SIMULATED  EXPLOSIONS  WITH  SENSITIVITY  TEST 


Problem 

Number 

Depth 

(m) 

Ms 

mr 

AT 

sec 

00 

Y— 1 1 

=»  o 

r-H 

X 

H0 

x IQ"18 

M2 

x IQ'18 

1* 

c 159 

4.38 

4.66 

★ 

0.1 

0.424 

■1 

1* 

c 159 

4.52 

4.66 

★ 

1.0 

0.304 

■S3 

0.611 

1(7 

c 159 

4.40 

4.66 

1.0 

0.304 

ESS 

0.558 

U 

c 159 

4.38 

4.66 

0.1 

0.424 

0.468 

0.631 

la 

4.59 

4.69 

★ 

0.1 

0.475 

0.378 

0.607 

lb 

c 207 

4.62 

4.69 

★ 

1.0 

0.524 

0.464 

0.700 

2 a 

4.60 

4.69 

1.0 

0.524 

0.558 

0.765 

3 a 

c 253 

4.29 

4.70 

★ 

0.1 

0.308 

0.378 

0.488 

3b 

c 253 

4.34 

4.70 

1.0 

0.371 

0.413 

0.555 

3a 

c 253 

4.38 

4.70 

1.0 

0.371 

0.451 

0.584 

10 

c 207 

4.11 

4.62 

★ 

0.1 

0.142 

0.108 

0.178 

12 

531 

4.17 

4.67 

★ 

0.1 

0.057 

0.206 

0.213 

26 

1000 

4.08 

4.70 

0.1 

0.057 

0.207 

0.215 

110 

207 

3.72 

4.69 

★ 

0.1 

0.052 

0.115 

0.126 

112 

531 

3.82 

4.70 

★ 

0.1 

0.048 

0.146 

0.154 

126 

1000 

3.86 

4.69 

fr 

0.1 

0.052 

0.174 

0.182 

c-- indicates  a crater  was  formed. 

★ 

--indicates  was  set  to  zero  on  last  time  step. 
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The  surface-wave  magnitude  falls  in  the  proper  range  for  150-KT 
explosion  events  when  compared  with  a list  of  measured  M$  versus  yield 
values  from  47  NTS  events  [11].  Of  the  47  events,  4 had  yields  between 
155  KT  and  170  KT,  with  values  from  3.97  to  4.32.  Our  results  for 
granite  are  on  the  high  side  of  that  distribution,  but  it  is  not  known 
whether  a granite  shot  is  included  in  the  NTS  list. 

The  PSR  M$  values  for  the  granite  shots  as  shown  in  Fig.  3 are  all 
higher  than  the  wet-sandstone  values,  but  it  is  felt  that  this  distinction 
does  not  have  conclusive  theoretical  support,  because  is  calculated 
from  force  alone  and  not  from  both  forces  and  displacements.  Richter's 
[12]  magnitude  expression  based  on  total  transmitted  energy--which  is 
proportional  to  the  log  of  the  energy  in  ergs--is  calculated  for  each  of 
the  explosion  events  and  presented  as  MR  in  Table  2.  Because  the  energy 
output  is  about  the  same  for  both  granite  and  sandstone  media--but  since 
granite  has  a larger  modulus  of  rigidity--one  would  expect  larger  forces 
for  granite.  But  because  the  rate  of  energy  involves  the  product  of 
force  and  rate  of  displacement  (velocity)  on  the  magic  circle,  one  would 
expect  larger  displacements  for  sandstone.  The  model  assumptions  used 
are  not  able  to  preserve  the  above-mentioned  distinction  between  the  rela- 
tive magnitudes  of  the  force  components  for  granite  and  sandstone.  This 
distinction  would  require  that  a pure  surface  wave  passing  through  the 
magic  circle  should  conserve  its  energy,  and  therefore  if  granite  forces 
were  larger  than  sandstone  forces  on  one  side  of  the  magic  circle,  the 
granite  forces  should  also  be  larger  on  the  other  side.  This  require- 
ment is  not  fulfilled,  however,  for  the  Green's  function  method  used  here, 
because  the  response  is  measured  in  centimeters  of  displacement  per  dyne 
of  force,  which  means  that  the  larger  forces  of  granite  will  produce  dis- 
placements that  are  larger  than  the  sandstone  values  of  displacement, 
contradicting  the  equivalent  conclusion  based  on  conservation  of  energy 
applied  to  the  two  media. 

The  above  ambiguity  of  the  model  does  not  apply  to  the  relative 
values  of  different  shots  in  the  same  medium.  For  granite,  the  variation 
of  is  caused  by  changes  in  the  shapes  of  the  time-dependent  forces 
acting  on  the  magic  circle,  and  not  changes  in  medium  parameters.  The 
slight  change  in  MR  as  a function  of  depth  indicates  that  some  mechanism 
is  withholding  more  energy  at  shallower  depths,  which  probably  represents 
energy  required  to  form  the  crater  and  propel  ejecta.  The  variation  of 
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M$  values  for  granite  suggests  significant  interaction  with  the  free  sur- 
face, and  also  suggests  that  this  interaction  has  a strong  enough  depth 
dependence  to  cause  a measurable  (0.3  units)  change  in  M$. 

Some  confirmation  is  given  to  the  values  for  granite  by  the  sensi- 
tivity test,  whose  range  of  values  is  shown  by  the  error  bars  in  Fig.  3, 
and  which  reveals  that  the  overall  shape  of  the  granite  curve  is  unchanged 
if  M$  is  assumed  to  lie  within  the  error  bars. 

The  interpolated  wet-sandstone  curve  of  Fig.  3 shows  no  strong  depth 
effect,  but  the  resolution  of  the  data  points  in  the  region  of  cratering 
(one  point)  is  not  enough  to  draw  any  conclusion.  The  two  deep  noncrater- 
ing wet-sandstone  shots  have  nearly  the  same  value,  which  is  to  be 
expected  according  to  simple  elastic  scaling  as  long  as  the  depth  is  much 
less  than  one  Rayleigh  wavelength. 

The  "elastic"  wet-sandstone  Ms  values  are  all  consistent  with  a reduced 
interaction  with  the  free  surface,  although  elastic  here  refers  only  to 
the  absence  of  cracking  and  gravity.  Generally,  the  higher  values  of 
for  the  nonelastic  wet  sandstone  seem  to  be  due  to  residual  block  motion 
(a  free-surface  interaction)  that  prolongs  the  force-time  histories  and 
allows  more  energy  into  the  20-sec  period  range.  If  one  follows  this  line 
of  reasoning,  the  shallowest  elastic  wet-sandstone  shot  would  have  the 
smallest  M$  value  because  the  significant  interaction  with  the  free  surface 
would  be  over  with  in  a shorter  time  period  and  therefore  would  deposit 
less  energy  in  long-period  components. 

Next,  all  of  the  results  presented  in  Table  2 and  Fig.  3 are  based  on 
time  histories  of  force  that  were  modified  in  a somewhat  ad-hoc  manner. 

The  original  force-time  histories  were  characterized  by  an  initial  sharp 
pulse,  followed  by  a rebound,  some  reflection  arrivals,  and  then  a tail 
of  low-intensity  oscillating  force.  The  cut-off  point  on  the  oscillating 
tail  would,  in  general,  be  at  a time  when  the  net  vertical  impulse  was  not 
zero;  and  this  small  residual  vertical  impulse  would  dominate  the  spectral 
amplitude  near  the  20-sec  periods.  Therefore,  after  about  2 sec,  at  the 
end  of  simulation,  one  extra  time  step  was  added,  with  a resultant  force 
that,  when  applied  for  AT  sec,  would  cause  the  net  vertical  impulse  to  be 
zero.  This  extra  force  at  the  end  is  not  completely  specified  by  the 
requirement,  VQ  = 0,  where 


21 


V0-/  FV(t)  dt  * (5) 

0 

The  unspecified  parts  are  the  distribution  of  the  vertical  force  over  the 

magic  circle  and  the  prescription  for  specifying  the  horizontal  component. 

The  distribution  of  the  vertical  component  is  taken  to  be  uniform  and  the 

horizontal  component  is  either  set  to  zero  everywhere  or  given  the  value 

of  the  previous  time  step.  This  choice  is  indicated  on  Table  2 (p.  18)  by 

an  asterisk  (*)  affixed  to  the  value  of  AT  when  the  horizontal  component, 

Fh,  is  set  to  zero,  and  by  no  indication  when  the  previous  values  of 

are  used.  Table  2 shows  various  combinations  of  aT  and  Fu  tried  for 

★ “ 

granite.  Only  the  combination  AT  = 0.1  is  used  for  wet  sandstone.  In 

either  instance,  the  velocities  for  the  last  time  step  are  set  to  zero, 

and  so  the  total  energy  transmitted  is  not  affected. 

Of  the  various  force  values  used  on  the  artificial  last  step,  all  M$ 

values  fall  within  the  error  bars  shown  on  the  granite  curve  of  Fig.  3. 

The  largest  difference  is  for  Problem  1,  the  shallowest  shot,  where  M$ 

varied  from  4.38  to  4.52,  a difference  of  0.14  units.  Figure  4 shows  some 

★ 

detailed  time  histories  for  Problem  1 for  the  case  AT  = 0.1  . Since  the 
net  vertical  force  (Fig.  4a)  must  integrate  to  zero,  the  spectral  source 
amplitudes  have  a lowest-order  dependence  on  Vj,  where 

T 

Vj(T)  = j Fv(t)  • t dt  . (6) 

o' 

One  can  clearly  see  the  artificial  last  0.1  sec  of  vertical  force  emerge 
above  the  oscillating  tail  (Fig.  4a)  to  drive  VQ  to  zero  at  the  end  time. 

A more  realistic  force-time  history  would  have  involved  two  oppositely 
directed  vertical-force  components,  one  upward,  due  to  buoyancy  while  the 
ejecta  is  in  free  fall;  the  other  downward,  due  to  small  impact.  The  final 
value  of  would  therefore  depend  on  momentum  carried  by  the  spall  and 
the  periods  of  the  different  possible  trajectories.  The  approximations 


M Oft?S 


e 


Fig.  4--Details  of  shallowest  granite 

simulation 

a.  Total  vertical  force  versus  time 

b.  Cumulative  vertical  Impulse 
versus  time 

c.  Cumulative  first-time  moment 
of  vertical  force  versus  time 

d.  Weighted  horizontal  force 
versus  time 

e.  Cumulative  horizontal  impulse 
versus  time 
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used  here  are  still  of  value,  however,  because  they  show  that  various 
reasonable  ways  to  terminate  affect  M$  less  than  the  depth  of  burial 
(DOB)  physical  parameter. 

The  magnitude  of  the  transient  impulse  due  to  ejecta  was  estimated 
to  be  about  one-tenth  the  value  due  to  spall  for  a typical  cratering 
I explosion  [13],  which  gives  a spall  inpulse  of  about  1018  dyne-sec  and 

an  ejecta  impulse  of  10^  dyne-sec.  These  transients  can  be  clearly 
seen  in  Fig.  4b,  except  some  interpretation  of  the  curve  must  be  made 
because  we  see  in  the  figure  the  reaction  of  the  impulse  on  the  rest  of 

• the  world,  and  there  seems  to  be  some  rebounding  taking  place.  The  value 

of  Vn  at  T = 2.1  sec,  just  before  the  last  artificial  time  step,  is 
u 1 7 

about  -10  dyne-sec,  which  can  be  interpreted  to  be  the  downward  impulse 
delivered  to  the  world  as  the  ejecta  is  thrown  up.  Additionally,  the 

• transient  spall  impulse  is  assumed  to  have  died  out. 

With  the  above  interpretations  of  the  shape  of  Fig.  4b,  we  find  good 
agreement  between  the  estimated  values  for  spall  and  ejecta  impulse,  and 
the  numerical  values  on  the  figure. 

• From  this  point  on,  since  the  ejecta  is  in  free  fall,  there  should 
be  a net  upward  force  due  to  the  lightening  of  the  load  removed  by  the 
free-falling  ejecta,  and  this  upward  force  would  cause  Vq  to  increase 
positively  in  value,  passing  through  zero  and  continuing  upward  until  the 

• ejecta  falls  back  to  earth,  at  which  time  Vq  returns  to  a value  of  zero. 
Since  this  imagined  time  history  was  not  included  in  the  source  data,  we 
do  not  know  for  certain  how  V^(t)  (Fig.  4c)  would  be  affected;  but  it  is 
not  expected  to  be  greatly  affected  because  of  the  results  of  the  sensi- 
tivity analysis  described  above. 

Finally,  some  results  of  a theoretical  nature  were  found  concerning 
the  correlation  of  V-j  and  Hq  with  M&,  where 

T 

H0  = / F*(t)  dt  ; (7) 

0 

that  is,  H-.  is  the  zero-th  time  moment  of  the  weighted  horizontal  force  F 
★ u 

Fu  Is  not  the  total  horizontal  force  on  the  magic  circle,  because  by 
n * 

symmetry  that  force  must  be  zero  at  each  instant  of  time.  Instead,  F^  is 
a measure  of  the  surface-wave-producing  ability  of  oppositely  directed 
horizontal  forces  separated  by  a distance  2r.  According  to  the  source 
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calculations  of  the  PSR  M$  model,  the  resultant  F^(t)  is  a summation 
over  the  individual  ring  sources  as  given  by 

FH(t)  = 2 i FHj(t)  * (8> 

j R 

where  is  the  x-coordinate  of  the  j-th  ring  force,  and  CR  is  the 

Rayleigh-wave  phase  velocity  for  20-sec  period  waves.  One  can  easily  > 

see  that  HQ  and  V1  both  have  units  of  (dyne-sec  ). 

In  the  full-blown  propagation  calculations,  each  Fu.(t)  and  each 

nj 

F\/ j ( t ) enters  into  the  seismometer  response,  with  factors  representing 
frequency  and  depth  dependencies;  but  most  of  this  dependence  is  only 
slowing  changing,  and  can  be  factored  out  to  a certain  degree  of 
approximation. 

After  this  approximation  is  made,  the  result  is  that  the  M value 

s 

is  proportional  to  M,,  where  the  subscript  2 indicates  that  it  is  the 
2 4 

coefficient  of  w in  the  spectral  expansion  of  the  teleseismic  ground 
response,  and  where 

M2  = Vvl  + Ho  • O) 

Tor  the  vertical  moment  , it  is  clear  that  the  net  vertical  force  Fy(t) 
must  reach  zero  after  some  time  period--say,  after  all  ejecta  have  re- 
turned to  earth.  The  formation  of  a cavity  at  the  explosion  site, 

★ 

however,  has  the  effect  of  leaving  a residual  outward  force  FR(t) 
acting  at  the  magic  circle,  and  so  HQ ( t ) will  continue  to  increase 
indefinitely  with  time,  as  can  be  seen  in  Figs.  4d  and  4e.  But  since 
the  time  moments  of  the  forces  are  used  only  as  coefficients  in  the 
spectral  expansion,  we  could  always  go  back  and  replace  the  expansion 
by  exp(-iwt),  which  would  serve  to  cut  off  the  artificial  growth  of 
the  individual  spectral  components.  For  the  present  effort,  however, 
the  expansion  must  suffice. 

Figure  5 shows  a scatter  plot  of  versus  for  the  entries  in 
Table  2 (p.  18).  Apparently,  a quite  good  correlation  exists  between 

> 


__ a. 


-Scatter  plot  of  M versus  M~  for  all  entries  In  Table  2 (p.  18)  showlnq 
a good  correlation.  c 
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the  single-source  parameter  and  the  resulting  value  of  M$,  which 
indicates  that  the  calculational  procedure  is  somewhat  transparent  to 
the  details  of  diffraction  from  the  composite  ring  source.  The  regres- 
sion line  of  Fig.  5 intersects  the  M axis  at  M = 3.9,  which  repre- 

. s s p 

sents  the  contribution  due  to  terms  of  higher  order  than  w in  the 
spectral  expansion.  The  three  "elastic"  wet-sandstone  shots  fall  in 
a well-defined  group  below  the  regression  line  of  the  nonelastic  shots. 
For  these  elastic  shots,  the  higher  order  terms  of  the  spectral  expansion 
must  be  lower  in  magnitude,  possibly  because  the  lack  of  cracking  stops 
the  late-time  block  motion  of  the  corresponding  nonelastic  shots  and 
therefore  less  energy  gets  deposited  in  the  long-period  waves. 

The  graphs  for  the  other  eight  shots  are  similar  in  form  to  those 
for  the  shallowest  granite  shot,  and  similar  interpretations  concern- 
ing ejecta,  etc.,  can  be  made  on  them.  Since  many  of  these  interpreta- 
tions would  be  redundant,  no  further  comment  is  made  on  them  here  (but 
they  all  are  presented  in  Appendix  G). 

CONCLUSIONS 

The  main  conclusion  of  this  section  is  that  the  surface  interaction 
with  an  underground  explosion  can  cause  a measurable  change  (0.3  units) 
in  the  teleseismic  surface-wave  magnitude,  M$,  for  different  depths  of 
burial,  at  least  for  the  simulated  explosions  in  granite.  The  sand- 
stone shots  were  not  done  in  the  appropriate  depth  ranges  to  allow  such 
a conclusion. 

However,  assuming  that  cracking  and  ejection  of  the  near-surface 
material  are  the  dominant  forms  of  free-surface  interaction  for  the 
wet  sandstone  shots,  which  is  supported  by  comparing  elastic  with  non- 
elastic results,  then  we  can  conclude  that  cracking  and  ejection  are 
the  dominant  source  variables  that  influence  the  shape  of  the  versus 
DOB  curves  for  both  granite  and  wet  sandstone  in  Fig.  3. 

We  have  also  found  that  the  surface-wave  Green's  function  approach 
to  the  surface-wave  propagation  calculation  can  be  made  internally  con- 
sistent and  that  it  is  a physically  justifiable  procedure. 
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III.  BODY  WAVES 


We  have  reached  three  main  conclusions  regarding  the  influence 
of  burial  depth  on  body  waves: 

1.  Particularly  as  it  influences  the  contribution  of  reflected 
waves  to  irrotational  motion,  burial  depth  must  be  taken  into 
account  in  relating  values  of  mb  to  yield.  In  wet  sandstone, 
m^  was  found  to  vary  by  0.38  (and,  hence,  estimated  yield  by 
a factor  of  2.6  to  3.0)  over  the  burial-depth-interval  from 
200m  to  550m. 

2.  For  depths  >200m,  magnitudes  that  are  almost  depth-independent 
can  be  obtained  from  the  "b"  phase  of  seismometer  displacement; 
the  role  of  such  magnitudes  (denoted  "m  ")  in  yield  estimation 

d 

is  at  present  minor,  and  should  be  enlarged. 

3.  If  reflection  is  assumed  to  take  place  elastically  at  the 
ground  surface,  then  estimates  of  the  effects  of  burial  depth 
on  body  waves  emerge  that  are  at  most  qualitatively  correct. 

To  compute  the  range  of  variation  of  mb,  or  the  depth-interval 
over  which  values  of  m^  span  that  range,  we  must  take  account 
of  both  gravity  and  tensile  cracking  on  the  motion  of  the 
medium. 

Qualitatively,  the  first  two  of  these  conclusions  are  not  new  (e.g., 
see  Ref.  14  and  AFTAC-TR-77-37) . Previously,  however,  they  have  been 
reached  by  assuming  reflection  to  be  elastic--a  procedure  that  we  find 
Inadequate  for  quantitative  purposes  (conclusion  3). 

BODY-WAVE  MAGNITUDES  VERSUS  SHOT  DEPTH:  INTERACTION  OF  DIRECT 
ANb  RELATED  waves 

Conclusions  1,  2,  and  3 above  follow  from  Fig.  6.  This  figure 
contains  both  our  computed  m^-values  for  inelastic  reflection,  and 
computed  values  that  show  how  m^  would  vary  if  tensile  cracking  and 
gravity  were  absent  (semielastic  reflection).  Curve  3 in  the  figure 
refers  to  purely  elastic  reflection,  treated  approximately.  As  is 
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evident  from  Fig.  6,  mb  does  not  vary  monotonical ly  with  burial  depth 
» in  any  of  the  three  cases.  Points  obtained  from  phase  b of  our  syn- 

thetic seismograms  (ma-values)  also  appear  in  Fig.  6;  their  relative 
independence  of  depth  is  clear. 

The  magnitudes  of  Fig.  6 are  determined  by  the  combined  effects 
i of  two  waves:  1)  a nearly  spherical  (and  initially  compressive)  wave 

driven  away  from  the  shot  point  by  high  pressures  that  attend  explosive 
energy  release;  2)  a reflected  wave  (mainly  expansive)  due  to  relief 
of  spherical-wave  stresses  at  the  ground  surface.  More  precisely,  for 
a 150-KT  explosion  at  a given  depth,  the  teleseismic  body-wave-field 
manifested  in  Fig.  6 can  be  viewed  as  the  sum  of  1 ) a "direct"  wave, 
defined  as  the  wavetrain  the  burst  would  generate  in  a whole-space; 
and  2)  a "relief"  wave  created  by  reflection  of  the  direct  wave  from 
the  ground  surface.  In  deriving  curve  3,  the  initial  wave  is  taken 
as  the  field  calculated  for  a burst  in  a uniform  whole-space  of  wet 
sandstone,  and  the  relief  wave  as  its  negative  (implosive)  image;  the 
centers  of  the  two  wave-fields  are  separated  by  twice  the  depth  of 
burial.  In  every  other  instance,  both  waves  are  generated  as  part  of 
a complete  numerical  solution  of  the  equations  of  motion. 

For  explosions  at  the  ground  surface,  the  direct  wave  and  relief 
wave  are  produced  at  almost  the  same  time,  and  hence  largely  cancel 
each  other.  As  burial  depth  increases  from  zero,  however,  the  relief 
wave  begins  to  lag  behind  the  direct  wave  (relief  starts  only  after 
the  direct  wave  reaches  the  ground  surface),  and  the  cancelling  effects 
of  one  wave  on  the  other  diminish.  At  a still  greater  burial  depth, 
the  relief  wave  arrives  just  as  the  negative  stage  of  direct-wave 
motion  (i.e.,  motion  back  toward  the  shot  point)  ends.  As  a result, 
at  points  on  nearly  vertical  downgoing  rays,  the  passage  of  time  is 
marked  by  a large,  continuous  change  from  downward  to  upward  displace- 
ment; when  thus  in  phase,  the  two  waves  reinforce  each  other.  Ulti- 
mately, shot  depths  are  reached  when  the  relief  wave  lags  so  far  be- 
hind the  direct  wave  that  the  two  no  longer  interact. 

The  effects  noted  have  counterparts  in  seismometer  response.  Thus, 
as  shot  depth  increases,  mfa  increases  at  first  and  then  decreases, 
arriving  finally  at  a whole-space  value  (^1250m  for  wet  sandstone). 


Before  reaching  its  whole-space  value,  m also  varies  with  burial  depth. 

a 

Since  ma  is  determined  by  relatively  early  features  of  seismograms, 
however,  the  relief  wave  ceases  to  have  any  effect  on  m at  a depth  of 

d 

burial  (^820  for  wet  sandstone)  and  a lag-time  that  are  relatively 
small.  For  reasons  noted  below  (p.  31),  relief-wave  effects  on  m are 
much  smaller  than  for  mb>  from  820m  to  depths  as  shallow  as  200m. 

FURTHER  COMMENTS  ON  FIG.  6:  ELASTIC  VERSUS  INELASTIC  REFLECTION 

Separating  the  far  field  of  a burst  into  a whole-space  ("direct") 
wave  and  a relief  wave  is  convenient  in  that  complete  body  waves  re- 
duce to  direct  waves  at  great  burial  depths.  More  important,  for 
bursts  of  a given  yield,  the  separation  in  question  assures  that  changes 
in  burial  depth  will  affect  the  relief  wave  alone.  In  addition,  given 
elastic  reflection,  the  shape  and  amplitude  of  the  relief  wave  itself 
become  independent  of  burial  depth  in  the  far  field.  The  teleseismic 
body  wave  then  depends  on  burial  depth  in  one  way  only:  the  time  by 
which  the  relief  wave  lags  behind  the  direct  wave  is  proportional  to 
that  depth. 

More  specifically,  seismograms  are  affected  by  lag-time.  Consider 
first  elastic  near-field  reflection.  As  shown  teleseismically,  dominant 
periods  in  the  direct  displacement  wave  run  from  2 sec  to  3 sec  (com- 
ponents with  shorter  periods  are  damped  relatively  strongly  by  the 
earth  as  they  run  from  near-field  to  far-field).  On  the  other  hand, 
the  AFTAC  seismometer  is  a narrow-band  instrument  most  sensitive  to 
displacement  harmonics  with  periods  of  'v-0.5  sec.  The  first  forced 
oscillation  of  the  seismometer  should  therefore  take  place  with  a 
period  somewhat  longer  than  0.5  sec.  Hence,  as  the  relief  wave's 
lag-time  varies  from  zero  to  a value  somewhat  larger  than  0.25  sec, 
the  signal  should  vary  from  complete  interference  with  the  direct  wave 
to  maximum  reinforcement.  As  a result,  given  elastic  reflection  at  a 
longitudinal  wavespeed  that  averages  2620m/sec  (from  shot  point  to 
ground-surface  and  back),  full  reinforcement  should  occur  at  a burial 
depth  somewhat  greater  than  330m.  As  it  happens,  our  synthetic  seismo- 
grams for  elastic  reflection  show  maximum  reinforcement  at  a depth  of 
about  450  m.  However,  factors  other  than  the  amplitude  of  seismometer 
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response  enter  the  definition  of  mb  (see  p.  35),  which  therefore  reaches 
I its  maximum  value  (assuming  elastic  reflection)  at  a depth  of  870m. 

For  a granite  medium  with  twice  the  longitudinal  wavespeed  of  wet  sand- 
stone, maximum  reinforcement  of  seismometer  response  (still  assuming 
elastic  reflection)  will  occur  at  a burial  depth  very  nearly  twice  450m, 
and  mb  will  probably  attain  its  maximum  at  2 * 870m  (or  1740m). 

If  we  take  into  account  cracking  and  the  earth's  gravitational  field, 
both  the  shape  and  amplitude  of  the  reflected  wave  at  teleseismic  dis- 
tances change  with  depth  of  burial.  The  head  of  the  relief  wave  still 
moves  downward  very  nearly  at  seismic  (elastic)  speed;  in  fact,  unload- 
ing proceeds  almost  elastically  until  cracking  occurs.  For  depths  of 
burial  down  to  1000m,  however,  the  period  of  elastic  unloading  is 
negligibly  short.  It  ends  when  the  medium  cracks,  a process  that  occurs 
computationally  as  soon  as  the  mean  stress  reaches  zero  (isotropic 
cracking).  Thus,  cracking  prevents  the  development  of  tension,  effec- 
tively limiting  the  amplitude  of  stress-relief  waves  to  overburden 
values.  Subsequent  recompression  of  the  medium,  due  to  its  reassembly 
under  gravity,  further  reduces  the  amplitude  of  the  relief  wave--at 
least  as  seen  along  rays  that  are  almost  vertically  downgoing  from  the 
shot  point.  Moreover,  since  the  motion  of  cracked  material  is  con- 
trolled largely  by  gravity,  development  of  the  relief  wave  is  slowed 
by  cracking:  instead  of  progressing  through  continuous  material  at 
elastic  wavespeed,  signals  are  repeatedly  interrupted  by  cracks. 

By  increasing  the  time  between  direct-  and  relief-wave  motions  at 
a given  depth,  gravity  and  cracking  have  the  general  effect  of  reducing 
the  shot-depth  interval  on  which  all  reinforcement-interference  phenomena 
occur.  Thus,  body-wave  magnitude  varies  more  rapidly  with  depth  in 
Fig.  6 when  gravity  and  cracking  are  taken  into  account  than  if  they  are 
not.  So  do  other  features  of  seismometer  response  (see  p.  35). 

The  combined  effects  of  cracking  and  gravity  are  greatest  for 
shallow  bursts.  The  farther  the  shot  point  from  the  ground  surface 
(increasing  depth  of  burial),  the  weaker  the  initial  wave  when  it 
arrives  at  the  surface  and  the  more  nearly  elastic  its  reflection  be- 
comes. Thus,  at  "great  enough"  shot  depths,  the  relief  wave  turns 
into  the  pP-wave  of  elastic-half-space  theory.  Table  3,  however,  shows 
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AMPLITUDE  AND  DURATION  OF  THE  RELIEF  WAVE  DUE  TO  A 150-KT 
BURST  IN  A WET-SANDSTONE  HALF-SPACE  VERSUS  DEPTH  OF  BURIAI 
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quantitatively  that  reflection  is  far  from  elastic  at  shot  depths  of 
200m  to  750m.  The  table  also  shows  that  reflection  becomes  increasingly 
elastic  as  depth  i ncreases--but  note  that  apart  from  the  smallest  of 
the  three  depths  tabulated,  higher  multipoles  of  the  field  probably 
affect  the  durations  and  amplitudes  listed.  Large  effects  of  gravity 
and  cracking  on  ref 1 ected-wave  motion  are  nevertheless  demonstrated  by 
Table  3,  since  the  "Elastic"  and  "Inelastic"  calculations  differ  in  no 
other  way  than  exclusion  or  inclusion  of  gravity  and  cracking.  Those 
same  effects  are  seen  in  still  more  detail  in  the  displacement  pulses 
of  Fig.  7.  In  particular,  note  that  at  the  shallowest  burial  depth, 
for  which  the  displacements  are  due  almost  entirely  to  monopole  radia- 
tion, hardly  a jot  remains  of  an  otherwise-prominent  pP-wave,  when 
cracking  and  gravity  are  taken  into  account. 

For  all  burial  depths  in  Fig.  7,  relatively  early  arrival  of  the 
first  disturbance  in  the  instance  of  gravity  and  cracking  is  due  to  an 
increase  in  p-wave  speed  with  increasing  overburden.  By  contrast,  in 
all  calculations  without  gravity  and  cracking,  the  sandstone  half-space 
is  considered  homogeneous;  its  uniform  p-wave  speed  is  set  to  the  value 
associated  with  overburden  at  a depth  of  531m.  Dispersion  in  a medium 
whose  wavespeed  increases  with  depth  also  accounts,  in  the  main,  for 
the  lower  peak  displacement  seen  at  each  depth  in  Fig.  7 for  calcula- 
tions with  gravity  and  cracking.  Resistance  to  cavity  growth  that 
increases  with  shot  depth  is  another  consequence  of  variable  overburden-- 
but  a minor  one--since  the  medium's  shear-strength  does  not  depend  on 
mean  stress.  For  a shot-depth  of  207m,  however,  shear  failure  was  ob- 
served in  the  relief  wave  as  well  as  the  direct  wave.  The  field  calcu- 
lated for  that  shot  depth  assuming  perfect  reflection  (see  Fig.  7 for 
pulses)  therefore  differs  from  the  other  two  fields  (i.e.,  with  gravity 
and  cracking,  and  without).  For  all  depths,  moreover,  perfect  reflec- 
tion eliminates  mode  conversion;  all  explosively  generated  energy  then 
resides  in  an  irrotational  field,  whence  relief-wave  amplitudes  are 
larger  for  perfect  reflection  than  if  we  merely  suppress  gravity  and 
cracking . 

Note  that  use  of  an  isotropic  model  of  cracking  in  the  calculations 
tends  to  minimize  crack  formation,  and  hence  to  reduce  associated 
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SHOT  DEPTH=207m 
RANGE=5290m 


SHOT  DEPTH=531m 
RANGE=4979m 


3 


SHOT  DEPTH=744m 
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Fig.  7--Displacement  pulses  in  the  elastic  near-field  of  a 150-kiloton 
burst,  for  three  depths  of  burial  in  a half-space  of  wet  sand- 
stone. The  pulses  refer  to  points  of  a down-going  ray  15°  off 
the  vertical,  at  distances  from  the  shot  point  denoted  by  values 
of  "RANGE." 
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gravitational  effects.  We  would  therefore  expect  greater  body-wave 
variations  than  those  of  Fig.  6,  Table  3,  and  Fig.  7 to  issue  from 
calculations  in  which  directional  cracking  is  permitted.  Also,  the 
kind  of  minor  seismogram  wiggle  that  leads  to  a discontinuous  change 
in  m^  in  Fig.  6 (p.  28)  would  almost  certainly  be  ignored  in  practice. 
Instead,  m^-values  would  be  obtained  from  a curve  like  1* — for  which 
our  calculated  variation  of  mb  with  depth  would  be  larger  than  0.38 
(the  value  quoted  above,  p.  27). 

DISCONTINUITIES  IN  m, 

b 

To  better  understand  the  results  shown  in  Fig.  6 (p.  28),  the 
specific  definition  of  mb  (used  in  preparing  the  figure)  must  be  con- 
sulted. In  brief,  of  the  several  terms  summed  to  give  m^,  only  one 
varies  when  depth  of  burial  is  changed.  That  term  is  log^g(A/TI), 
where,  as  calculated  for  Fig.  6,  A is  the  difference  between  the  most 
negative  value  of  a seismometer's  displacement;  the  following  maximum 
value,  T,  is  the  time  between  the  extremes  of  displacement  represented 
in  A;  and  I is  an  instrument  correction  that  depends  on  T (and  only  T) . 
Alternatively,  A could  be  interpreted  as  the  maximum  swing  of  a seismo- 
meter over  all  pairs  of  consecutive  extremes  of  its  displacement. 
Fortunately,  if  m^  is  defined  in  that  way,  no  changes  occur  in  curves 
1,  1’,  and  2 of  Fig.  6,  and  no  important  changes  in  curve  3.  However, 
if  A is  taken  as  the  difference  between  the  most  negative  seismometer 
displacement  and  the  most  positive  subsequent  displacement,  then  curve 
1 continues  along  1'  instead  of  dropping  discontinuously. 

A term  like  log,g(A/TI)  appears  in  almost  all  definitions  of  body- 
wave  magnitude  (though  often  with  I h 1),  and  leads  directly  to  the 
sharp  break  seen  (Fig.  6)  in  the  m^-depth  curve  for  inelastic  reflec- 
tion at  a depth  of  ^575m,  and  in  both  curves  for  elastic  reflection  at 
^865m.  Those  breaks  are  true  discontinuities.  They  occur  because  A, 

T,  and  (hence)  I,  are  discontinuous  functions  of  depth.  In  turn,  jumps 
appear  in  A and  T (hence,  I)  as  functions  of  burial  depth  because  the 
consecutive  extremes  of  seismometer  displacement  that  determine  A and 
T can  and  do  change  abruptly  with  depth.  In  their  turn,  abrupt  changes 


occur  in  those  extremes  because  the  total  response  of  the  seismometer 

(a  linear  operator  on  ground  displacement)  is  the  sum  of  its  individual 

responses  to  the  direct  wave,  and  to  the  relief  wave,  generated  by  any 

near-surface  burst  (p.  29).  Depending  on  the  delay  between  the  two 

waves--and,  hence,  on  burial  depth--that  sum  may  show  all  the  extremes 

of  each  individual  response;  alternatively,  a peak  or  trough  of  one 

response  curve  can  overlie  a peak  or  trough  of  the  other,  so  that  the 

two  are  fused  into  a single  extreme  of  the  complete  seismogram.  Thus, 

the  number  of  extremes  itself  depends  on  burial  depth.  As  a result, 

it  is  almost  inevitable  that  sudden  changes  appear  in  the  amplitude,  A, 

of  any  given  cycle  of  seismometer  oscillation;  for  example,  a "second 

maximum"  whose  value  is  used  to  compute  A at  one  burial  depth  may  be 

the  third  maximum  at  another  depth.  The  approximate  method  by  which 

curve  3 of  Fig.  6 was  constructed  will  now  serve  to  show  more  concretely 

how  discontinuities  in  m.  arise. 

b 

SPECIFICS  OF  INTERFERENCE  AND  REINFORCEMENT:  APPROXIMATION 

Of  P£EEect  reflection 

Increasing  the  depth  of  a burst  of  given  yield  causes  the  relief 
wave  originating  at  the  ground  surface  to  trail  farther  and  farther  be- 
hind the  direct  wave  always  created  in  such  events  (p.  29).  In  fact, 
as  required  by  causality,  the  effects  of  reflection  are  delayed,  rela- 
tive to  initial  seismometer  deflection,  by  a time  equal  to  that  of  the 
first  reverberation  between  shot  point  and  ground-surface.  Those 
effects  are  reviewed  in  detail  here  assuming  reflection  to  be  perfect. 
The  relief  wave  then  becomes  the  negative  image,  across  the  ground 
plane,  of  the  direct  wave;  also,  the  ground  plane  becomes  a surface 
dividing  a uniform  whole-space  into  two  half-spaces. 

Take  first  the  instance  of  great  lag  (very  deep  burial),  in  which 
direct-wave  displacements  determine  m^.  As  a function  of  time  at  fixed 
position,  those  displacements  form  a pulse  that  looks  like  a positive 
half-cycle  of  a sine  wave  (positive  displacement  being  away  from  the 
shot-point);  a minor  oscillation,  due  to  spherical  divergence  of  the 
field,  also  appears  as  the  pulse  ends.  Accordingly,  the  seismometer 
at  first  moves  away  from  the  point  of  burst.  Being  an  oscillator  with 
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its  own  characteristic  period  M sec),  it  thereafter  exhibits  stages 
of  both  positive  and  negative  displacement  as  it  swings  about  its  equi- 
librium position.  Being  a damped  oscillator,  it  soon  comes  nearly  to 
rest.  Thus,  for  a long  while,  the  seismogram  contains  only  two  posi- 
tive lobes,  and  one  negative  lobe,  of  consequence.  Eventually,  however, 
the  negative  of  that  response  appears,  as  the  seismometer  moves  under 
the  influence  of  the  relief  wave  reflected  frjm  the  ground  surface. 

The  individual  seismometer  responses,  and  their  superoosition  to  give 
its  total  response,  are  shown  in  Fig.  3,  along  with  both  near-  and  far- 
field  pulses  that  define  the  displacement  of  ground  material.  For  very 
deep  burial,  the  curve  of  total  response--seismogram  "A"  of  Fig.  8--is 
duplicated  when  the  individual  responses  are  plotted.  Simply  put,  the 
two  waves  do  not  interact  either  as  input  to  the  seismometer  or  as 
transformed  by  it. 

A like  set  of  plots  is  also  shown  in  Fig.  8 for  a shot  depth  of 
890m.  At  that  depth,  individual  responses  to  the  direct  wave  and  to 
the  relief  wave  do  interact.  The  initial  response  of  the  seismometer 
to  the  relief  wave,  which  is  negative,  occurs  at  such  time  as  to  inter- 
fere with  its  second  positive  outward  swing  in  response  to  the  direct 
wave.  A minor  oscillation  is  thereby  introduced  into  the  latter  re- 
sponse curve.  At  a slightly  shallower  burial  depth  (860m),  and  hence 
smaller  lag-time,  overlap  of  the  initial  negative  relief-wave  response 
with  the  second  positive  direct-wave  response  is  reduced.  As  a result, 
the  low-amplitude  maximum-minimum  oscillation  of  seismogram  3 is  de- 
graded to  a dimple  in  a curve  of  constantly  increasing  displacement 
(seismogram  C of  Fig.  8). 

The  change  from  seismogram  B to  seismogram  C has  evidently  taken 
place  continuously  with  depth.  It  is  equally  clear,  however,  that  the 
difference  between  the  all-time  minimum  displacement  and  next  maximum 
(A.j,  say)  has  changed  discontinuously,  and  with  it  the  time  between 
those  two  extremes.  So  has  the  maximum  displacement-swing  found  among 
all  pairs  of  consecutive  extremes  (A^,  say);  that  swing  occurs  between 
points  5 and  6 of  seismogram  B (point  5,  the  earlier  of  the  two,  is  a 
maximum),  and  points  2 and  3 of  seismogram  C (the  earlier  point,  2,  is 
a minimum).  Thus,  if  A1  or  A^  is  taken  as  the  amplitude  A,  and  T is 


NEAR  FIELD  FAR  FIELD  TWO  RESPONSES  TOTAL  RESPONSE 


Fig.  8--Ground  motion  and  seismometer  response:  150-kiloton  burst  in  a 

wet-sandstone  half-space;  perfect  ground-surface  reflection;  three 
burial  depths.  From  left  to  right,  the  four  curves  for  each  shot 
depth  show,  vs.  time,  1)  displacement  of  material  along  a down- 
going  ray  15°  off  the  vertical,  100  kilometers  (km)  from  the  shot, 
2)  displacement  on  that  ray  at  5000  km  ("FAR  FIELD"),  3)  separate 
far-field  seismometer  responses  to  the  direct  wave  and  relief  wave 
("TWO  RESPONSES"),  and  4)  total  seismometer  response. 
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equated  to  the  time  between  the  extremes  that  determine  A,  then  the 
term  log-j Q( A/TI ) will  jump  on  varying  the  shot  depth  from  890m  (case  B) 
to  860m  (case  C).  In  deriving  curves  1,  2,  and  3 of  Fig.  6 (p.  28), 
a definition  of  A^type  was  chosen  for  A (i.e.,  from  all-time  minimum 
displacement  to  the  next  maximum),  which  led  to  the  discontinuities  seen 
in  the  figure.  On  the  other  hand,  switching  to  a definition  of  A^-type 
would  leave  curves  1 and  2 unchanged,  as  well  as  the  depth  at  which 
curve  3 jumps.  As  just  noted,  however,  the  size  of  the  jump  would 
change  in  the  instance  of  perfect  reflection  (curve  3 of  Fig.  6);  points 
2 and  3 of  seismogram  B determine  m^  if  A is  A-j-like,  and  points  5 and 
6 if  A is  A2~ 1 i ke . 

A transition  from  the  pattern  of  seismometer  deflection  seen  in 
seismogram  C to  that  of  seismogram  B is  evident  in  Fig.  9.  Figure  9 
contains  plots  like  those  of  Fig.  8 ("Two  Responses"  excepted),  but 
applies  to  our  full-blown  calculations  of  motion  with  gravity  and  crack- 
ing. For  the  corresponding  calculations  without  gravity  and  cracking, 
only  pattern  C appears  (Fig.  10);  at  burial  depths  less  than  ^870m,  some 
mechanism  must  act  to  enhance  lag-time,  as  gravity  and  cracking  do 
(p.  33),  or  else  direct-  and  relief-wave  responses  will  not  be  separated 
enough  to  produce  pattern  B. 

As  burial  depth  is  reduced  from  860m  in  the  instance  of  perfect 
reflection,  the  initial  negative  relief-wave  response  (R^,  say)  over- 
laps more  and  more  fully  the  first  negative  direct-wave  response  (D^). 
Since  has  more  than  twice  the  amplitude  of  R^,  however,  only  a 
modest  increase  (30  to  40  percent)  takes  place  in  the  swing  between  the 
first  minimum  of  seismometer  deflection  (point  2 of  seismograms  A,  B, 

C),  and  either  the  next-earlier  or  next-later  maximum  (point  1 or  point 
3).  Most  of  the  latter  amplitude  increase  with  decreasing  depth  (i.e., 
between  points  2 and  3)  is  due  to  reinforcement  of  the  second  positive 
direct-wave  response  by  the  first  positive  relief-wave  response. 

Major  changes  occur  in  the  value  of  A that  enters  mb  only  when 
burial  depths  become  so  shallow  (<400m)  that  the  first  positive  relief- 
wave  response  starts  to  cancel  its  image  (the  first  negative  direct- 
wave  response).  As  shot  depth  is  further  decreased,  the  amplitude  A 
drops  at  an  ever-increasing  rate  and,  with  it,  m^  itself.  The  same 


Fig.  9--Near-  and  far-field  displacement  pulses,  and  seismometer  responses,  for 
150-kiloton  bursts  at  three  depths  in  a wet-sandstone  half-space,  calcu- 
lated uith  gravity  and  cracking.  The  pulses  refer  to  points  of  a down- 
going ray  15°  off  the  vertical,  at  near-field  distances  from  shot  point 
denoted  by  "RANGE";  the  far  field  is  5000-km  away. 
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Fig.  10--Near-  and  far-field  displacement  pulses,  and  seismometer  responses, 
for  150-kiloton  bursts  at  three  depths  in  a wet-sandstone  half- 
space, calculated  without  gravity  and  cracking.  The  pulses  refer 
to  points  of  a down-going  ray  15°  off  the  vertical,  at  near-field 
distances  from  shot  point  denoted  by  "RANGE";  the  far  field  is 
5000-km  away. 
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process  of  interference  also  causes  m and  its  amplitude-factor  to  fall. 

a 

That  decay,  however,  occurs  at  smaller  shot  depths  (<200m)  than  for  mb: 
it  takes  more  complete  overlap  of  the  relief  wave  and  direct  wave  to 
bring  about  decay  of  the  first  minimum  of  total  response  than  to  effect 
decay  of  the  second  maximum;  to  increase  the  waves'  overlap,  shot  depth 
must  be  reduced. 

THE  FACTORS  THAT  ENTER  m,  AND  m.  , AND  THEIR  VARIATION  WITH  DEPTH 

a D 

Figure  8,  and  the  discussion  just  concluded  of  reinforcement  and 

interference  of  direct  and  relief  waves,  serves  to  explain  how  the 

amplitude-factor  A varies  with  depth,  whether  for  m or  m.  . The  varia- 

a b 

tion  of  m^  (or  mg)  with  depth  (Fig.  6)  is  not  determined  by  A alone, 
however,  but  by  log^(A/TI)  (p.  35).  To  make  that  variation  intel- 
ligible, A,  T,  I,  and  A/TI  have  been  plotted  in  Fig.  11  (in  semilog 
fashion)  versus  depth  of  burial.  Values  of  A,  etc.,  are  separately 
displayed  in  Fig.  11  for  four  instances.  In  two  instances,  gravity 
and  cracking  are  included  in  our  calculations  of  motion;  the  values 

that  enter  m are  shown,  and  so  are  the  values  that  enter  m,  . Then, 
a D 

for  the  same  calculations,  but  without  gravity  and  cracking,  the  fac- 
tors A,  etc.,  are  plotted;  in  addition,  the  figure  contains  mb-factors 
for  the  instance  of  perfect  reflection  from  the  ground  surface. 

Figure  11  shows  that  the  factors  A,  etc.,  that  enter  rn  vary  little 

d 

with  depth.  The  same  should  therefore  be  true  of  m , and  it  is  (Fig.  6, 

a 

p.  28).  Also,  as  explained  above  (p.  35),  A is  a weak  function  of  depth 

for  m.  as  well  as  m , except  at  quite  shallow  depths  (<200m).  On  the 
d a 

other  hand,  the  period  T that  enters  mb  does  change  appreciably  with 
depth  of  burial.  However,  the  most  striking  result  conveyed  by  the 
figure  is  that  for  mb>  the  instrument  correction,  I,  varies  a good  deal 
more  with  depth  than  either  A or  T.  Thus,  with  some  justification,  the 
variation  of  mb  with  depth  in  Fig.  6 might  be  viewed  as  a figment  of 
the  mb*definition.  In  particular,  as  discussed  below,  most  of  the 
depth-dependence  of  mb  would  disappear  if  the  term  1 og-j gA/T I were  re- 
placed by  1 og i qA  in  defining  m . Yet,  that  view  is  not  exhaustive: 
the  dependence  of  T on  shot  depth,  from  which  the  variation  of  I follows, 
is  real  (though  as  already  noted,  and  as  Fig.  11  shows,  gravity  and 
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Fiq.  11 --Three  depth-dependent  factors  that  determine  body-wave  magnitudes  ma 
and  m^,  versus  shot  depth,  for  150-kiloton  bursts  in  a wet-sandstone 
half-space.  Amplitude  A is  equal  to  the  seismometer  deflection  from 
its  all-time  minimum,  either  to  its  next-earlier  maximum  (ma),  or  to 
its  next-later  maximum  (mh).  T is  twice  the  time-interval  Between 
the  extremes  that  define  A;  the  instrument  correction  I varies  with 
T only.  Changes  in  magnitude  with  depth  are  identical  to  changes  in 
log1n(A/TI)--also  shown. 
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cracking  must  be  taken  into  account  to  calculate  it).  The  strong 
dependence  of  I on  depth,  and  more  directly  the  changes  in  T that  it 
signifies,  therefore  have  a positive  side.  Those  T-changes,  which  are 
amplified  by  gravity  and  cracking,  may  offer  a means  of  deducing  shot 
depth  from  the  early  portions  of  seismograms,  whereas  the  same  data 
provide  values  of  1 og -j q A from  which  to  determine  explosive  yield--a 
subject  to  which  we  now  turn. 

CORRELATION  OF  BODY-WAVE  SOURCE  STRENGTH  WITH  FEATURES  OF  SEISMOGRAMS: 
FURTHER  COMMENTS  ~ ~~~  " 

Discontinuities  like  those  of  Fig.  6 (p.  28)  are  troubling  chiefly 
because  they  mean  that  m^  has  been  so  defined  as  to  measure  different 
things  at  different  depths.  Although  can  easily  be  redefined  to 
avoid  jumps,  it  is  no  simple  matter  to  obtain  an  internally  consistent, 
useful  measure  of  body-wave  strength  in  the  process.  For  example,  A 
could  be  taken  as  the  swing  from  all-time  minimum  seismometer  displace- 
ment to  the  greatest  subsequent  maximum,  as  in  the  case  of  points  2 and 
5 of  seismogram  B (Fig.  8,  p.  38).  Indeed,  in  a real  seismogram,  an 
oscillation  like  that  of  points  3 and  4 in  seismogram  B would  almost 
certainly  be  ignored,  and  the  A-factor  reckoned  instead  from  points  2 
and  5.  At  greater  depths,  however,  points  2 and  5 of  seismogram  B 
become  2 and  5 of  A.  Thus,  in  an  actual  seismogram,  5 would  be  mixed 
with  late-arriving  signals  of  all  kinds,  making  the  true  body-wave 
amplitude  at  that  point,  and  even  the  point's  identification,  highly 
uncertain.  Moreover,  the  amplitude  of  the  oscillation  marked  by  points 
3 and  4 would  then  have  become  large,  and  points  2 and  3 would  again 
span  the  full  amplitude-range  covered  by  the  seismogram  in  its  early 
cycles  of  pure  body-wave  response. 

Alternatively,  normalized  weights  w1  and  w 2,  respectively,  might 
be  used  to  combine  linearly  two  different  seismometer  deflections: 

1)  the  swing  S-|  from  an  all-time  minimum  deflection  m1  to  the  next 
maximum  M^ ; 2)  the  swing  S.,  from  m^  to  the  greatest  subsequent  maximum. 
Weight  w^  could  be  made  a smooth  function  of  the  time  T^  taken  to  com- 
plete swing  S2,  with  the  value  of  the  function  equal  to  zero  when  T£ 
exceeds  a specified  multiple  of  the  period  covered  by  S-| ; w-j  would 
vanish  smoothly  as  zero  is  approached  by  the  time  from  maximum  M1  to 
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I the  next  minimum  (e.g.,  from  point  3 to  point  4 of  the  seismogram  B 

in  Fig.  8).  As  we  see  it,  little  of  value  is  likely  to  emerge  from 

such  algebraic  thrashing  about,  because  it  does  nothing  to  solve  the 

, basic  problem  of  untangling  the  direct  wave  and  relief  from  one  another. 

L.  For  example,  if  a weighted  average  like  that  described  above  were  used 

to  compute  the  factor  A,  then  m,  would  still  measure  two  different 

b 

things  at  different  depths--al bei t with  a smooth  transition  from  one 
thing  to  the  other. 

On  balance,  given  its  near-freedom  from  relief-wave  influence  for 
all  but  cratering  bursts,  m^  seems  more  closely  related  than  m^  to  body- 
wave  strength--yet  it  too  has  drawbacks.  For  one,  the  amplitude  of  the 
first  significant  positive  swing  of  the  seismometer  is  low  relative  to 
those  of  its  subsequent  swings,  and  is  often  lost  in  background  noise. 

In  most  instances,  the  factor  A is  therefore  appreciably  smaller  for 
ma  than  m^,  though  not  by  more  than  a factor  of  two.  As  a result,  the 
minimum  yield  for  detection  of  explosions  is  somewhat  higher  when  m 

3 

is  used  rather  than  m^.  More  serious  are  the  large  effects  of  back- 
, ground  noise  on  the  period  T that  enters  m„  , as  opposed  to  m.  --effects 

' magnified  by  the  instrument  correction,  I.  For  example,  noise  that 

makes  the  first  maximum  uncertain  by  a factor  of  two  may  only  alter  A 
by  10  percent  to  20  percent  while  giving  rise  to  a 40-percent  error 
in  T and  a still  larger  error  in  the  product  TI.  Moreover,  the  first 
maximum  of  seismometer  response  is  likely  to  be  more  sensitive  to 
variations  in  site  medium  than  are  its  later  extremes  of  deflection. 

The  difficulties  cited  lead  us  to  suggest  that  the  A-factor  alone 
be  used  in  place  of  A/TI  as  a measure  of  body-wave  strength.  Indeed, 
according  to  Figs.  6 and  11  (pp.  28  and  43),  relief-wave  effects  on 
m^,  including  its  discontinuity,  are  much  larger  than  on  the  value  of 
log  A that  enters  m^.  We  therefore  consider  log^A  a better  index  of 
body-wave  strength  than  the  m^  quantity  to  which  it  contributes.  A 
combination  of  that  value  of  log^A  with  ma  (or  the  A-factor  associated 
therewith)  might  be  more  useful  still. 
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